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Abstract 

We derive and analyze a simplified formulation of the numerical viscosity terms 
appearing in the expression of the numerical fluxes associated to several High- 
Resolution Shock-Capturing schemes. After some algebraic pre-processing, we give 
explicit expressions for the numerical viscosity terms of two of the most widely used 
flux formulae, which implementation saves computational time in multidimensional 
simulations of relativistic flows. Additionally, such treatment explicitely cancells 
and factorizes a number of terms helping to amortiguate the growing of round-off 
errors. We have checked the performance of our formulation running a 3D relativis- 
tic hydrodynamical code to solve a standard test-bed problem and found that the 
improvement in efficiency is of high practical interest in numerical simulations of 
relativistic flows in Astrophysics. 
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The numerical study of the evolution of multidimensional relativistic flows turns out to 
be a topic of crucial interest in, at least, two different scientific fields: Nuclear Physics 
(studies of the properties of the equation of state for nuclear matter via comparison 
of simulations and experiments of heavy ion collisions) and Relativistic Astrophysics. 
The field of Numerical Relativistic Astrophysics is recently undergoing an extraordinary 
developement after the important efforts of people working in building up robust codes 
able to describe many different astrophysical scenarios, such that relativistic jets in quasars 
and microquasars, accretion onto compact objects, collision of compact objects, stellar 
core collapse and recent models of Gamma- Ray bursts (see, e.g., the recent review in 
[7] and references therein). Thus, the improvement in the efficiency of multidimensional 
hydro-codes becomes a necessity. 

It is well known the performance of modern high-resolution shock- capturing techniques 
(HRSC) in simulations of complex classical flows. Most of the HRSC methods are based 
on the solution of local Riemann problems (i.e., initial value problems with discontiu- 
ous initial data) and since 1991 [12] several Riemann Solvers or Flux Formulae have been 
specifically designed in relativistic fluid dynamics (see, e.g., [10,7] for a review on Riemann 
solvers in Relativistic Astrophysics). In addition, in a recent paper [13] we showed the 
way for applying special relativistic Riemann solvers in General Relativistic Hydrodynam- 
ics, hence any future new Riemann solver, exhaustively analyzed in Special Relativistic 
Hydrodynamics (SRH), could be applied to get the numerical solution of local Riemann 
problems in General Relativistic Hydrodynamics. Consequently, the interest of the results 
we obtain in this note goes beyond the domain of SRH and can be easily extended to 
General Relativistic Hydrodynamics. 

For consistency, we start by summarizing the basics of the HRSC techniques. A system 
of conservation laws [8] is 



where uG 9ft d is the vector of unknowns and f l (u) is the flux in the i-direction. In the 
above system (1) we can define a 5 x 5-Jacobian matrix B l (u) associated to the flux in 
the i-direction as: 



The system is said to be hyperbolic if the Jacobian matrices have real eigenvalues. 
The main ingredients of a HRSC algorithm are: 

i) A finite discretization of the equations in conservation form (1). Using a method of 
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lines, this discretization reads: 
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where subscripts i, j, k are related, respectively, with x, y and ^-discretizations, and refer to 
cell-centered quantities. The cell width, in the three coordinate directions are, respectively, 
Ax, Ay and Az. 

ii) Quantities f i+ i j k , gj J+ i fc and h ijifc+ i are called the numerical fluxes at the cell inter- 
faces. These numerical fluxes are, in general, functions of the states of the system at each 
side of the cell interface. Some HRSC methods derive expressions for the numerical fluxes 
by giving a consistent flux formulae or solving local Riemann problems, with an exact 
[11] or approximate Riemann solver, after a cell reconstruction procedure that gives the 
state at both sides of the interface, denoted by L (left state) and R (right state). Several 
monotonic cell reconstruction prescriptions have been given in the scientific literature to 
achieve different orders of spatial accuracy [15,16,9]. 

For clarity, from now on we will omit the indexes relative to the grid and restrict our 
study to the xi-splitting of the above system (1), assuming that the vector of unknowns 
satisfies u = n{x\,t). 

We have focussed our analysis to some of the most popular HRSC algorithms, and ana- 
lyzed their expressions for the numerical fluxes. Hence, the sample considered is: HLLE 
[6,3], Roe [14], Marquina (M) [2], and a modified Marquina's flux formula (MM) [1]. The 
above selection gathers the most fundamental differences among the large sample of HRSC 
flux formulae. HLLE is the simplest one, it does not need the full spectral decomposition 
of the Jacobian matrices. Roe's solver linearizes the information contained in the spectral 
decomposition into an average state. Marquina's (and its modified version) flux formula 
considers the information coming from each side of a given interface (it is not a Riemann 
solver) and, in some astrophysical applications, has produced better results in modelling 
complex flows. 

After some algebraic work, all these flux formulae can be cast into the following general 
form: 

f{u L ,u R ) = l - ((J + J L )f L + (J- l R ){ R + (Q L u L - Q R u R )) (4) 

where f L,R stands for the fluxes evaluated at the states u L ' R and J is the unit matrix. Fol- 
lowing Harten [5] , the Q c ' n terms in the above equation will be referred as the numerical 
viscosity matrix. 

Matrices I L,R and Q L > R can be expressed as linear combinations of the projectors onto 
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TABLE I 

Parameters in the numerical fluxes. 



Flux b p c p 

Formulae 
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Table 1 

In the above table we have introduced the quantities = max(0, \ R , \+) and = 
min(0, X R , A_), A + and A_ are, respectively, the maximum and minimum of X p , a p = 

max (j A^ |, | Xp \j and j3 p = \ (sgn{\p) + sgn(Xp)j. We denote by u the state of the system 
according to Roe's average. 

each eigenspace, i.e., the direct product of the corresponding left and right eigenvectors 
l p ,r p associated to the p-th characteristic field (p=l,...,d), 



P =i 
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where superscripts L, R indicate that the eigenvectors are evaluated at the state u L ' R . 
The values of the coefficients b p and c p appearing in the above definitions of matrices I L > R 
and Q L ' R depend on the eigenvalues A p as shown in Table I, for the four flux formulae 
analyzed. 

Several comments concerning Table I are in order: 

i) If we take into account the orthonormality relations between the right and left eigen- 
vectors 

X>p = 1 (7) 

P =i 
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and the fact that the coefficients b p and c p are, in the case of HLLE, independents of p, 
then matrices X L,R and Q L ' R are, trivially, the unit matrix multiplied by the corresponding 
factors. 

ii) For HLLE's and Roe's flux formulae their corresponding matrices X L ' R and Q L ' R satisfy 
the relations: T L = 1 R , Q L = Q R 

iii) As it is well known, the knowledge of the spectral decomposition of the Jacobian 
matrices is a basic ingredient to build up Riemann solvers or many flux formulae. Never- 
theless, while HLLE's flux formula only needs the values of the maximum and minimum 
speeds of propagation of the signals, Roe's and Marquina's flux formulae need explicitly 
the full knowledge of the spectral decomposition, including right and left eigenvectors. 

The system governing the evolution of multidimensional relativistic perfect fluids can be 
written in Cartesian coordinates in the form (1), with d — 5, where, in units such that 
the speed of light c = 1, the vector of unknowns u is given by 

u=(D,S\ S 2 , S\ r) T , (8) 



the fluxes are defined by 

P = (Dv\ SV + p8 u , S 2 v l + P 5 2i , SV + P S 3i , S i - Dv'Y (9) 



where D(= pW) is the rest mass density, S^(= phW 2 v^) is the j-component of the mo- 
mentum density, and r(= phW 2 — p — pW) is the energy density, W — (1 — v 2 ) -1 / 2 is 
the Lorentz factor, p is the rest-mass density, p the pressure and h the specific enthalpy 
given byh = l + e + p/p with e being the specific internal energy. The system of SRH is 
closed with an equation of state p = p(p,e) from which the local sound speed, c s , can be 
obtained 

^ = | + ^)|. no) 



In a previous paper [2] we derived the explicit analytical expressions for the full (right 
and left) spectral decomposition. We denote the five characteristic fields by p — 1, . . . , 5 = 
— , 0, 0, 0, +, in the standard ordenation. 

A very worthy simplification on the calculation of matrices Q arises when some eigenvalue 
is degenerate, i.e., when the system is not strictly hyperbolic. In SRH, like in multidimen- 
sional Newtonian hydrodynamics, there is a linearly degenerate field, p — 0, such that the 
corresponding eigenvalue Ao is triple (the system in the j-direction splitting is not strictly 
hyperbolic, although the set of eigenvectors is complete). According to equation (6), and 
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using the ort honor mality relations between the right and left eigenvectors 

ErSyk = *™-r?j;-r™J» (11) 
k=i 

where m, n = 1, . . .5 denote the components of the 5-vector, it is possible to eliminate 
the three eigenvectors associated to the degenerate field and to write down the following 
simplified expression (omitting L,R superscripts) 

Q mn = ^mn + ^ + _ C() ) r ™/« + ( c _ _ Co ) r ™P. (12) 

Notice that only r± and l± are needed to evaluate the numerical viscosity. The same pro- 
cedure can be applied to any system of conservation laws where one of the eigenvectors has 
multiple degeneracy, because orthogonality relations always allow us to skip the explicit 
dependence on one of the vector subspaces of the spectral decomposition. In particular, 
it could be of great interest in the case of the equations of relativistic magnetohydro- 
dynamics where, in the ansatz of a directional splitting, similar degeneracy arises in the 
structure of the characteristic fields associated to each one of the fluxes. 

The explicit formulae for the numerical viscosity term corresponding to the system of 
equations of special relativistic hydrodynamics are: 

HLLE's flux formulae. Since the numerical viscosity matrix is proportional to the identity, 
the application of the above recipes is obvious. 

Roe 's flux formulae. The numerical flux across some given interface can be written 

f(u L ,u*) = i[f L + f R + q] (13) 

q being the five-vector calculated from the corresponding numerical viscosity matrices of 
Table I: 

q = Q( U L - u R ) = QAu (14) 

In Roe's Riemann solver the quantities relative to the spectral decomposition are evaluated 
using the corresponding Roe-averages of the left and right states, denoted by u (see [14], for 
the Newtonian case and [4] for the relativistic case). In practice, other particular averaging 
(e.g., arithmetic means) have also been used. Note that in the following expressions (15) 
all quantities are evaluated using Roe's average, except for Au m . After some algebra, the 
viscosity vector associated to the numerical flux in the j-direction is 

qi = | A | Ami + Xa 
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q2 = 

93 = 

q 4 = 
1h = 



A | Au 2 + hW(v x Xa + XbSjx) 
A | Am 3 + hW(vyXa + Xb^jy) 
A | Aw 4 + hW(v z x a + Xb^'z) 

A | Am 5 + /lW(Xa + VjXb) - Xa 



where 
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(17) 
(18) 



M and MM- formulae. The numerical flux across a given interface can be written like 
equation (13) with 



q = q L - q fi 

^L,R _ qL,R^L,R 



(19) 
(20) 



Omitting the superscripts L, i? and taken into account the expressions in Table I for MM 
and the results in [2], the viscosity vector in the x-splitting is: 

h 2 

qt R = ^{M [A-CL+ - A+Q-] + p(c+K+ - c_N_)} + 
W f K v 2 + v 2 } 



c p- 



h \K-1 l-v 2 x 



q L, R = WW_ { MA+A _ [ fi+A+ _ fi _A_] + p{c + \ + A+K + - c_A_^_N_)} + 

q L, R = hW_^ {M [n+A _ _ n _ A+] + p(c+N+ _ c K )} + 
f W^ 2 l + 2^ 2 (^ 2 + t; 2 )l 

c n^T + } 

q L, R = h_W Vz {M _ n _ A+ ^ + p(c+K+ _ c K )} + 

W^ 2 l + 2W^ 2 + t; 2 )l 
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h 2 

q^ R = — {M - A+n_V_] + p[c + K + V + - c_K_T?_]} + 



W ( hW-JC (2W-l)(tj , .... , 

c oP^r v 1 + 1 — > ( 21 ) 




A 

' // [ a: i 1 i - v 

with the following auxiliary quantities 

M = phW 2 (JC-l), n± = c±(v x -\ T ), V± = hWA±-l, 



y. K l BP 

K — Cg p oe 



A± = 1 - ^ 



1 - v x A± 



A = /i 3 W/(/C - 1)(1 - u x u a .)(A h A+ - ^_A_) 
N± = ± {-«x - ^ V - ^)(2/C - 1)(^ - AfcA±) + /C^l ± A ± } (22) 

where quantities c± ; o are given in Table I and A is the determinant of the matrix of 
right-eigenvectors. 

The corresponding viscosity vectors in the other directions are trivially obtained by a 
cyclic permutation of subindices x, y, z. 

We have tested the efficiency of our numerical proposal, for Roe's and MM's flux formulae, 
running GENESIS (a 3D special relativistic hydro-code [1]), without any optimization at 
compilation level, in a SGI Origin 2000. A standard initial value problem has been chosen: 
Pl = 10, 6l = 2, vl — 0, 7l = 5/3, pr — 1, €r — 10~ 6 , vr — and 7^ = 5/3, where the 
subscript L (R) denotes the state to the left (right) of the initial discontinuity. This test 
problem has been considered by several authors in the past (see [1] for details in ID, 2D 
and 3D). 

We have compared the performance of two different implementations of the numerical 
flux subroutine: 

i) Case A, stands for the results obtained using our analytical prescription. This means 
to write down, in the numerical flux routine, just the expressions derived here for the 
viscosity vector q. 

ii) Case B, stands for the results obtained running the code with a standard high-efficiency 
subroutine for inverting matrices (we use a LU decomposition plus an implicit pivoting 
which is, for general matrices, 0(N 3 )). This subroutine is called to get the left eigenvectors 
from the matrix of right eigenvectors and is adapted to the particular dimensions of the 
matrices (3 x 3, in ID, 4 x 4, in 2D and 5 x 5, in 3D). Hence, unlike case A, now we 
have to calculate numerically the following quantities: the matrix of left eigenvectors, the 
characteristic variables and, finally, the components of the viscosity vector q. 

Table II summarizes the results: the direct implementation of our numerical viscosity 
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TABLE II 

CPU time (in microseconds). 









TCI (fis) 






Roe 


MM 


Case 


# Zones 


Case A 


Case B 


Case A Case B 


ID 


100 x 1 x 1 


12.2 


53.8 


23.8 118.9 


2D 


20 x 20 x 1 


25.5 


181.8 


49.0 373.5 


3D 


14 x 14 x 14 


39.4 


431.9 


75.7 879.0 



Table 2 

Time per numerical cell and iteration (TCI) in microseconds employed by the numerical flux 
routine in our test-bed problem, for three different grids. 

formulae leads to an improvement of the efficiency (in terms of CPU time) of the numerical 
fluxes subroutine in a factor which, in 3D calculations, ranges between about eleven and 
twelve depending on the particular flux formula used. When comparing Roe's and MM's 
cases a factor two -in favour of Roe- arises due to the fact that MM's flux formulae needs 
to compute two viscosity vectors (one per each side of a given interface), unlike Roe's 
flux formula which needs only one viscosity vector evaluated at the average state. As it 
must be, the efficiency increases with the number of spatial dimensions involved in the 
problem due to the computationally expensive matrix inverting operations performed at 
each interface to get the numerical fluxes. Since the numerical flux routine is, typically, 
one of the most time-consuming, it translates into a speed up factor between two and four 
in the total execution time, depending on the specific weight of the flux formulae routine 
in each particular application. 

Our formulation also gives a unified description of the numerical fluxes (4), permitting 
a unique implementation with the possibility of switching in cases when the utilisation 
of a specific flux formula is more appropriate. In addition, due to the fact that we have 
eliminated, in the generalized MM's flux formula, all the conditional clauses, the efficiency 
is ensured either for scalar or vectorial processors. 

Another worthy by-products of our algebraic pre-processing concerns with the significant 
reduction of round-off errors, as a consequence of the number of operations suppressed 
and factorization. One of the important issues in designing a multidimensional hydro- 
code is the accurate preservation of any symmetries present in a physical problem. A 
numerical violation of these symmetries could arise as a consequence of accumulation 
of round-off errors in the calculation of the numerical fluxes, as we have explained in a 
previous paper [1]. The algebraic simplifications, shown in the present paper, reduce the 
number of operations and cure such problem. 
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Two last additional consequences arise from our work. First is that similar expressions 
can be worked out for any non-linear hyperbolic system of conservation laws for which the 
full spectral decomposition is known. In particular, when some of the vectorial subspaces 
has multiple degeneracy, a similar algebraic preprocessing is very convenient. The other 
important consequence is that an appropriate combination of a simplified formulation of 
the numerical viscosity together with the use of special relativistic Riemann solvers in 
General Relativistic Hydrodynamics [13], should allow a very easy and efficient extension 
to General Relativistic Hydrodynamics. 
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